# generate heatmap and return the ordered matrix for other plot to follow the same order
#when input a list of matrix, it will just clustering based on the first matrix, then the others will also be print out as the first's order
# 
# Author: yaping
###############################################################################
library(gplots)
library(colorRamps)
library(RColorBrewer)
library(GenomicRanges)
##hclust in this package is much faster, and allow more rows...
library(fastcluster)
source("/home/uec-00/yapingli/code/mytools/R/heatmap.3.R")
#source("/Volumes/HD_2/Documents/workspace/dnaase_google_script/edu.usc.epigenome.dnaase.Rscript/heatmap.3.R")
#source("/Volumes/storage/hpcc/uec-00/yapingli/code/mytools/R/heatmap.3.R")

jet.colors <-colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan","#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
jet.rev.colors <-colorRampPalette(c("#7F0000","red","#FF7F00","yellow","#7FFF7F","cyan","#007FFF","blue","#00007F"))
white2red<-colorRampPalette(c("white","red"))
cyan2red<-colorRampPalette(c("cyan","grey","red"),bias=2)
white2green<-colorRampPalette(c("white","darkgreen"))
blue2yellow<-colorRampPalette(c("blue","yellow"))
subClusterColors<-c("red","orange","blue","black","cyan","purple","palegreen3","yellow","lightpink2","ivory4")
#if subCluster is 1, then average plot would only generate one single plot.
#otherwise, it will provide average plot in different sub clusters.

generateHeatmapForNOMeSeq<-function(fileName,prefix="heatmap", scale=1000, move_step=20,bin_size_align=1, rowOrder=NULL, pngOutput=FALSE,pdfOutput=TRUE, regionToCluster=c(38:62), dendgramAlone=FALSE, averagePlot=FALSE, 
		subClusterNum=1, subClusterOrder=NULL,inputClusterOrder=TRUE, printSideBar=F,y_min=0, y_max=100, outputSubClusterCordinate=F, outSideDendgram=NULL,maxNAallow=8,
		capLimit=FALSE, capUplimit=NULL, capDownLimit=NULL,logScale=FALSE, autoScale=FALSE, orderByMaxOccDisToTss=FALSE,MaxOccWindow=120,
		rowSideFiles=NULL, colSideFiles=NULL,maxGap=250, colToUse=rbind(c("cyan","white"),c("green","white"),c("orange","white"),c("red","white"),c("blue","white"),c("purple","white"),c("black","white")),heatMapCols=white2green(75), ...){
	content<-read.table(fileName,sep="\t",header=F)	
	content<-content[!duplicated(content),]
	#rownames(content)<-paste(content[,1],content[,2],content[,3],sep=":")
	rownames(content)<-paste(content[,1],content[,2],content[,3],content[,4],sep=":")
	#dataSeq<-seq(5+((length(content[1,])-5)/2)-scale, 5+((length(content[1,])-5)/2)-scale+scale * 2, by=move_step)
	dataSeq<-seq(5+((length(content[1,])-5)/2)-as.integer(scale/bin_size_align), 5+((length(content[1,])-5)/2)+as.integer(scale/bin_size_align), by=as.integer(move_step/bin_size_align))
	valueGch<-array()
	
	for(i in dataSeq){
		#valueGch<-cbind(valueGch,as.numeric(rowMeans(content[,(i-move_step/2):(i+move_step/2-1)], na.rm=T)))
		if(floor(i+move_step/(2*bin_size_align)-1)-ceiling(i-move_step/(2*bin_size_align)) <= 0){
			valueGch<-cbind(valueGch,as.numeric(content[,ceiling(i-move_step/(2*bin_size_align))]))	
		}else{
			#tmp<-content[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)]
			#x<-ifelse(rowSums(is.na(tmp))<(move_step-1),mean(tmp,na.rm=T),NA)
			#valueGch<-cbind(valueGch,x)
			valueGch<-cbind(valueGch,as.numeric(rowMeans(content[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)], na.rm=T)))	
		}
		
	}
	valueGch<-valueGch[,2:length(valueGch[1,])]	
	rownames(valueGch)<-rownames(content)
	
	if(autoScale){
		y_min=min(valueGch,na.rm =T)
		y_max=max(valueGch,na.rm =T)
		##use colMeans will be more suitable for autoscale
		
	}
	if(capLimit){
		if(is.null(capUplimit)){
			capUplimit=quantile(unlist(valueGch),probs=seq(0,1,0.05),na.rm =T)[20] #cap top 5%
			
		}
		if(is.null(capDownLimit)){
			capDownLimit=quantile(unlist(valueGch),probs=seq(0,1,0.05),na.rm =T)[2] #cap down 5%
			capDownLimit=ifelse(capDownLimit<0,0,capDownLimit) ##avoid some very small negative value close to 0...
		}
		y_min=capDownLimit
		y_max=capUplimit
		valueGch[valueGch<capDownLimit]=capDownLimit
		valueGch[valueGch>capUplimit]=capUplimit
		content[content<capDownLimit]=capDownLimit
		content[content>capUplimit]=capUplimit
		#valueGch<-apply(valueGch, 1:2, function(x) ifelse(x<capDownLimit,capDownLimit,ifelse(x>capUplimit,capUplimit,x)))
		#content[,5:length(content[1,])]<-apply(content[,5:length(content[1,])], 1:2, function(x) ifelse(x<capDownLimit,capDownLimit,ifelse(x>capUplimit,capUplimit,x)))

	}
	
	
	if(logScale){
		psedoCount=quantile(unlist(valueGch),probs=seq(0,1,0.01),na.rm =T)[2]
		y_min=log2(y_min+psedoCount) ##need to deal with y_min=0....
		y_max=log2(y_max+psedoCount)
		valueGch<-log2(valueGch+psedoCount)
	}
	
	#clustering step
	if(is.null(rowOrder)){
		if(orderByMaxOccDisToTss){ ##order by max Occupancy score positions' distance to TSS
			rowMaxOccScore=rep(2000,length(valueGch[,1]))
			rowMaxOccScorePos=rep(median(regionToCluster),length(valueGch[,1]))
			for(pos in regionToCluster){
				meanVector<-rowMeans(valueGch[,(pos-as.integer(MaxOccWindow/(move_step*2))):(pos+as.integer(MaxOccWindow/(move_step*2)))],na.rm=T)
				tmp<-meanVector-rowMaxOccScore
				rowMaxOccScore[!is.na(tmp) & tmp<0]=meanVector[!is.na(tmp) & tmp<0]
				rowMaxOccScorePos[!is.na(tmp) & tmp<0]=rep(pos,length(tmp[!is.na(tmp) & tmp<0]))
			}
			valueGch<-valueGch[order(rowMaxOccScorePos),]
			rowOrder<-rownames(valueGch)
			subClusterOrder=rep(1,length(valueGch[,1]))
			names(subClusterOrder)<-rowOrder
			
		}else{  ##do clustering rather than ordering
			#maxNAallow=quantile(rowSums(is.na(valueGch[,regionToCluster])),probs=seq(0,1,0.05),na.rm =T)[19]
			valueTmpGchNoCluster<-valueGch[rowSums(is.na(valueGch[,regionToCluster]))<=maxNAallow,]
			rownames(valueTmpGchNoCluster)<-rownames(valueGch[rowSums(is.na(valueGch[,regionToCluster]))<=maxNAallow,])
			#valueTmpGchNoCluster<-valueGch[rowSums(is.na(valueGch[,40:60]))<=5,]
			#rownames(valueTmpGchNoCluster)<-rownames(valueGch[rowSums(is.na(valueGch[,40:60]))<=5,])
			#ward, mcquitty,  average, 
			x<-hclust(dist(valueTmpGchNoCluster[,regionToCluster]),method="ward")
			valueGch<-valueTmpGchNoCluster
			rowOrder<-x$order
			
			if(is.null(subClusterOrder) & subClusterNum > 1){
				subClusterOrder<-cutree(x, k = subClusterNum)	
				inputClusterOrder=FALSE
			}
		}
		
	}else{
		rowOrder<-rowOrder[rowOrder %in% rownames(valueGch)]
	}
	
	##plot average plot, may be seperated into different sub-cluster order
	if(averagePlot){
		axisSeq<-seq(0-scale, scale, by=move_step)
		axisSeqForPlot<-seq(0-scale, scale, by=scale/4)
	}
	if(subClusterNum > 1 & averagePlot){
		pdf(paste("AveragePlot.",subClusterNum,"clusters.",fileName,".pdf",sep=""))
		for(clusterNum in c(1:subClusterNum)){
			cluster_name<-names(subClusterOrder[subClusterOrder==clusterNum])
			valueTmpCluster<-content[cluster_name,]
			valueGch1<-array()
			for(j in dataSeq){	
				#valueGch1<-cbind(valueGch1,mean(colMeans(valueTmpCluster[,(j-move_step/2):(j+move_step/2-1)], na.rm=T), na.rm=T))
				if(floor(j+move_step/(2*bin_size_align)-1)-ceiling(j-move_step/(2*bin_size_align)) <= 0){
					valueGch1<-cbind(valueGch1,mean(valueTmpCluster[,ceiling(j-move_step/(2*bin_size_align))], na.rm=T))	
				}else{
					#tmp<-valueTmpCluster[,ceiling(i-move_step/(2*bin_size_align)):floor(i+move_step/(2*bin_size_align)-1)]
					#x<-ifelse(colSums(is.na(tmp))<length(tmp[,1])-1,mean(tmp,na.rm=T),NA)
					#valueGch1<-cbind(valueGch1,mean(x,na.rm=T))
					valueGch1<-cbind(valueGch1,mean(colMeans(valueTmpCluster[,ceiling(j-move_step/(2*bin_size_align)):floor(j+move_step/(2*bin_size_align)-1)], na.rm=T), na.rm=T))	
				}
			}
			valueGch1<-valueGch1[,2:length(valueGch1[1,])]
			#valueGch1<-colMeans(valueTmpCluster,na.rm=T)			
			if(logScale){
				valueGch1<-log2(valueGch1)
			}
			
			
			plot(axisSeq,valueGch1,type="l",axes=FALSE,ylim=c(y_min,y_max),xlab="",ylab="",col=subClusterColors[clusterNum],lty=1,font=2,lwd=3)
			par(new=T)
		}
		axis(1,at=axisSeqForPlot,lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2)
		axis(2,at=seq(y_min,y_max,by=(y_max-y_min)/5),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2)
		title(paste("AveragePlot.",subClusterNum,"clusters.","\n",fileName,sep=""), cex.main = 0.6, font.main= 4, col.main= "black",xlab="Distance to center (bp)", ylab="Methylation")
		abline(v=0)
		dev.off()
	}
	
	valueTmpGch<-valueGch[rowOrder,]
	
	##generate row side colorbar
	colBar=NULL
	
	if(printSideBar){
		if(subClusterNum>1){
			subClusterColBar=subClusterColors[subClusterOrder]
			if(inputClusterOrder){
				names(subClusterColBar)=names(subClusterOrder)
			}
			subClusterColBar=subClusterColBar[rowOrder]
			colBar<-cbind(colBar,subClusterColBar)
		}
		if(!is.null(rowSideFiles)){
			for(fileOrder in c(1:length(rowSideFiles))){
				rowSideData<-read.table(rowSideFiles[fileOrder],sep="\t",header=F)
				rowSideData_loc<-GRanges(seqnames=rowSideData[,1],ranges=IRanges(rowSideData[,2],rowSideData[,3]),strand="*")
				print(rowSideFiles[fileOrder])	
				print(dim(rowSideData))
				valueTmp_loc_vector<-strsplit(rownames(valueTmpGch),":")
				valueTmpGch_loc<-GRanges(seqnames=sapply(valueTmp_loc_vector,function(x) x[1]),ranges=IRanges(sapply(valueTmp_loc_vector,function(x) as.integer(as.numeric(x[2])+as.numeric(x[3]))/2),width=1),strand="*")
				print(dim(valueTmp_loc_vector))
				print(dim(valueTmpGch))
				rowSideData_bar<-ifelse(countOverlaps(valueTmpGch_loc,rowSideData_loc,maxgap=maxGap)>0,colToUse[fileOrder,1],colToUse[fileOrder,2])
				colBar<-cbind(colBar,rowSideData_bar)
			}
			
		}else{
			colBar<-cbind(colBar,rep("white",length(subClusterColBar)))
			
		}
	}
	addAverageToHeatmap=FALSE
	if(averagePlot)
		addAverageToHeatmap=TRUE
	
	#dendgram
	if(dendgramAlone){
		fileNameDend=paste("Dendgram.",fileName,".pdf",sep="")
		pdf(fileNameDend)
		if(is.null(outSideDendgram)){
			plot(x)
			outSideDendgram<-as.dendrogram(x)
		}else{
			plot(outSideDendgram)
		}
		
		dev.off()
	}
	
	#output each cluster's coordinates
	if(outputSubClusterCordinate){
		for(clusterNum in c(1:subClusterNum)){
			cluster_name<-names(subClusterOrder[subClusterOrder==clusterNum])
			fileNameClust=paste("Cluster",clusterNum,".",fileName,".bed",sep="")
			clusterCor<-strsplit(cluster_name,":")
			clusterBed<-NULL
			for(num in c(1:length(clusterCor))){
				strand<-sub("NONE", ".", clusterCor[[num]][4])
				strand<-sub("NEGATIVE", "-", strand)
				strand<-sub("POSITIVE", "+", strand)
				cor<-c(clusterCor[[num]][1],clusterCor[[num]][2],clusterCor[[num]][3],".",".",strand)
				clusterBed<-rbind(clusterBed,cor)
			}
			
			write.table(clusterBed,fileNameClust,sep="\t",quote=F,col.names=F, row.names=F)
			
		}
	}
	
	#plotting step
	if(pdfOutput){
		fileNamePdf=paste(fileName,".pdf",sep="")
		pdf(fileNamePdf, width=5,height=7.5)
	}else{
		fileNamePng=paste(fileName,".png",sep="")
		png(fileNamePng)
	}

if((is.null(outSideDendgram) & !dendgramAlone)){
	rowV=FALSE
	rowD="none"
}else{
	rowV=outSideDendgram
	rowD="row"
}
#print(ifelse((is.null(outSideDendgram) & !dendgramAlone), FALSE, outSideDendgram))
		heatmap.3(valueTmpGch,
				Rowv=rowV,
				Colv=FALSE,
				dendrogram= rowD,
				distfun = dist,
				hclustfun = hclust,
				na.rm=TRUE,
				labRow = "",
				labCol = "",
				#labCol = seq(0-scale,scale,by=step),,
				cexCol = 0.05 + 1/log10(dim(valueTmpGch)[2]),
				#key=TRUE,
				keysize=1,
				trace="none",
				density.info=c("none"),
				margins=c(1, 1),
				RowSideColors=colBar,
				col=heatMapCols,
				#na.color=par("bg"),
			addAverage=addAverageToHeatmap,
			averageData=content,
			logScale=logScale,
			clusterNums=subClusterNum,
			clusterOrder=subClusterOrder,
			colAveragePlot=subClusterColors,
			xAxisSeqForAvePlot=axisSeqForPlot,
			axisSeq=axisSeq,
			dataStep=as.numeric(move_step),
			dataSeq=dataSeq,
			y_min=y_min,
			y_max=y_max,
			bin_size_in_align=as.numeric(bin_size_align),
				...
		)
		dev.off()

	capUplimit=y_max
	capDownLimit=y_min
	returnResult<-list(rowNames=rownames(valueTmpGch), subClusterOrder=subClusterOrder, capUplimit=capUplimit, capDownLimit=capDownLimit)
	return(returnResult)

	
}

##regionToCluster in all files input here are geather together to do clustering, which is good for identified out differentiate clusters out.
##only the rows exist in all rows will be kept, the row names are 1st file's coordinate.
generateHeatmapForNOMeSeqByDifferentiateClustering<-function(fileName,prefix="heatmap",sampleLen=101, scale=1000, step=20,bin_size_align=1, rowOrder=NULL, pngOutput=FALSE,pdfOutput=TRUE, fileNamesForClustering=NULL, regionToCluster=c(38:72), maxNAallow=8, dendgramAlone=FALSE, averagePlot=FALSE, 
		subClusterNum=1, subClusterOrder=NULL,inputClusterOrder=TRUE, printSideBar=F,yscale=100,
		rowSideFiles=NULL, colSideFiles=NULL,maxGap=250, colToUse=rbind(c("cyan","white"),c("green","white"),c("orange","white"),c("red","white"),c("blue","white"),c("purple","white"),c("black","white")),heatMapCols=white2green(75), ...){
	
	dataForClustering<-NULL
	commonRowsNames<-NULL
	if(length(fileNamesForClustering)>=2){
		for(num in c(1:length(fileNamesForClustering))){
			content<-read.table(fileNamesForClustering[num],sep="\t",header=F)
			rownames(content)<-paste(content[,1],content[,2],content[,3],sep=":")
			dataSeq<-seq(5+((length(content[1,])-5)/2)-as.integer(scale/bin_size_align), 5+((length(content[1,])-5)/2)+as.integer(scale/bin_size_align), by=as.integer(step/bin_size_align))
			valueGch<-NULL
			
			for(i in dataSeq){
				#valueGch<-cbind(valueGch,as.numeric(rowMeans(content[,(i-step/2):(i+step/2-1)], na.rm=T)))
				if(floor(i+step/(2*bin_size_align)-1)-ceiling(i-step/(2*bin_size_align)) <= 0){
					valueGch<-cbind(valueGch,as.numeric(content[,ceiling(i-step/(2*bin_size_align))]))	
				}else{
					valueGch<-cbind(valueGch,as.numeric(rowMeans(content[,ceiling(i-step/(2*bin_size_align)):floor(i+step/(2*bin_size_align)-1)], na.rm=T)))	
				}
				
			}
			#valueGch<-valueGch[,2:length(valueGch[1,])]
			rownames(valueGch)<-rownames(content)
			if(is.null(commonRowsNames)){
				commonRowsNames<-rownames(valueGch)
				
			}else{
				commonRowsNames<-commonRowsNames[commonRowsNames %in% rownames(valueGch)]
				valueGch<-valueGch[commonRowsNames,]
				dataForClustering<-dataForClustering[commonRowsNames,]
			}
			dataForClustering<-cbind(dataForClustering,valueGch[,regionToCluster])
			rownames(dataForClustering)<-commonRowsNames
		}
	}
	
	
	content<-read.table(fileName,sep="\t",header=F)	
	#rownames(content)<-paste(content[,1],content[,2],content[,3],sep=":")
	rownames(content)<-paste(content[,1],content[,2],content[,3],sep=":")
	#dataSeq<-seq(5+((length(content[1,])-5)/2)-scale, 5+((length(content[1,])-5)/2)-scale+scale * 2, by=step)
	dataSeq<-seq(5+((length(content[1,])-5)/2)-as.integer(scale/bin_size_align), 5+((length(content[1,])-5)/2)+as.integer(scale/bin_size_align), by=as.integer(step/bin_size_align))
	valueGch<-array()
	
	for(i in dataSeq){
		#valueGch<-cbind(valueGch,as.numeric(rowMeans(content[,(i-step/2):(i+step/2-1)], na.rm=T)))
		if(floor(i+step/(2*bin_size_align)-1)-ceiling(i-step/(2*bin_size_align)) <= 0){
			valueGch<-cbind(valueGch,as.numeric(content[,ceiling(i-step/(2*bin_size_align))]))	
		}else{
			valueGch<-cbind(valueGch,as.numeric(rowMeans(content[,ceiling(i-step/(2*bin_size_align)):floor(i+step/(2*bin_size_align)-1)], na.rm=T)))	
		}
		
	}
	valueGch<-valueGch[,2:length(valueGch[1,])]	
	rownames(valueGch)<-rownames(content)
	
	if(is.null(dataForClustering)){
		dataForClustering<-valueGch[,regionToCluster]
		commonRowsNames<-rownames(valueGch)
	}
	
	#clustering step
	if(is.null(rowOrder)){
		valueTmpGchNoCluster<-dataForClustering[rowSums(is.na(dataForClustering))<=maxNAallow,]
		rownames(valueTmpGchNoCluster)<-rownames(dataForClustering[rowSums(is.na(dataForClustering))<=maxNAallow,])
		#valueTmpGchNoCluster<-valueGch[rowSums(is.na(valueGch[,40:60]))<=5,]
		#rownames(valueTmpGchNoCluster)<-rownames(valueGch[rowSums(is.na(valueGch[,40:60]))<=5,])
		#ward, mcquitty,  average, 
		x<-hclust(dist(valueTmpGchNoCluster),method="ward")
		#valueGch<-valueTmpGchNoCluster
		rowOrder<-x$order
		
		if(is.null(subClusterOrder) & subClusterNum > 0){
			subClusterOrder<-cutree(x, k = subClusterNum)	
			inputClusterOrder=FALSE
		}
	}else{
		rowOrder<-rowOrder[rowOrder %in% commonRowsNames]
	}
	
	##plot average plot, may be seperated into different sub-cluster order
	if(subClusterNum > 0 & averagePlot){
		axisSeq<-seq(0-scale, scale, by=step)
		axisSeqForPlot<-seq(0-scale, scale, by=scale/4)
		pdf(paste("AveragePlot.",subClusterNum,"clusters.",prefix,".",fileName,".pdf",sep=""))
		for(clusterNum in c(1:subClusterNum)){
			cluster_name<-names(subClusterOrder[subClusterOrder==clusterNum])
			valueTmpCluster<-content[cluster_name,]
			valueGch1<-array()
			for(j in dataSeq){	
				#valueGch1<-cbind(valueGch1,mean(colMeans(valueTmpCluster[,(j-step/2):(j+step/2-1)], na.rm=T), na.rm=T))
				if(floor(j+step/(2*bin_size_align)-1)-ceiling(j-step/(2*bin_size_align)) <= 0){
					valueGch1<-cbind(valueGch1,mean(valueTmpCluster[,ceiling(j-step/(2*bin_size_align))], na.rm=T))	
				}else{
					valueGch1<-cbind(valueGch1,mean(colMeans(valueTmpCluster[,ceiling(j-step/(2*bin_size_align)):floor(j+step/(2*bin_size_align)-1)], na.rm=T), na.rm=T))	
				}
			}
			valueGch1<-valueGch1[,2:length(valueGch1[1,])]
			#valueGch1<-colMeans(valueTmpCluster,na.rm=T)			
			
			plot(axisSeq,valueGch1,type="l",axes=FALSE,ylim=c(0,yscale),xlab="",ylab="",col=subClusterColors[clusterNum],lty=1,font=2,lwd=3)
			par(new=T)
		}
		axis(1,at=axisSeqForPlot,lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2)
		axis(2,at=seq(0,yscale,by=yscale/5),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2)
		title(paste("AveragePlot.",subClusterNum,"clusters.",prefix,"\n",fileName,sep=""), cex.main = 0.6, font.main= 4, col.main= "black",xlab="Distance to center (bp)", ylab="Methylation")
		abline(v=0)
		dev.off()
	}
	
	valueTmpGch<-valueGch[rowOrder,]
	
	##generate row side colorbar
	colBar=NULL
	
	if(printSideBar){
		if(subClusterNum>1){
			subClusterColBar=subClusterColors[subClusterOrder]
			if(inputClusterOrder){
				names(subClusterColBar)=names(subClusterOrder)
			}
			subClusterColBar=subClusterColBar[rowOrder]
			colBar<-cbind(colBar,subClusterColBar)
		}
		if(!is.null(rowSideFiles)){
			for(fileOrder in c(1:length(rowSideFiles))){
				rowSideData<-read.table(rowSideFiles[fileOrder],sep="\t",header=F)
				rowSideData_loc<-GRanges(seqnames=rowSideData[,1],ranges=IRanges(rowSideData[,2],rowSideData[,3]),strand="*")
				
				valueTmp_loc_vector<-strsplit(rownames(valueTmpGch),":")
				valueTmpGch_loc<-GRanges(seqnames=sapply(valueTmp_loc_vector,function(x) x[1]),ranges=IRanges(sapply(valueTmp_loc_vector,function(x) as.integer(as.numeric(x[2])+as.numeric(x[3]))/2),width=1),strand="*")
				rowSideData_bar<-ifelse(countOverlaps(valueTmpGch_loc,rowSideData_loc,maxgap=maxGap)>0,colToUse[fileOrder,1],colToUse[fileOrder,2])
				colBar<-cbind(colBar,rowSideData_bar)
			}
			
		}
	}
	addAverageToHeatmap=FALSE
	if(averagePlot)
		addAverageToHeatmap=TRUE
	
	#plotting step
	if(pdfOutput){
		fileNamePdf=paste(prefix,".",fileName,".pdf",sep="")
		pdf(fileNamePdf, width=5,height=7.5)
		#pdf(fileNamePdf)
		heatmap.3(valueTmpGch,
				Rowv=FALSE,
				Colv=FALSE,
				dendrogram= c("none"),
				distfun = dist,
				hclustfun = hclust,
				na.rm=TRUE,
				labRow = "",
				labCol = "",
				#labCol = seq(0-scale,scale,by=step),,
				cexCol = 0.05 + 1/log10(dim(valueTmpGch)[2]),
				#key=TRUE,
				keysize=1,
				trace="none",
				density.info=c("none"),
				margins=c(1, 1),
				RowSideColors=colBar,
				col=heatMapCols,
				#na.color=par("bg"),
				addAverage=addAverageToHeatmap,
				averageData=content,
				clusterNums=subClusterNum,
				clusterOrder=subClusterOrder,
				colAveragePlot=subClusterColors,
				xAxisSeqForAvePlot=axisSeqForPlot,
				yAxisSeqForAvePlot=seq(0,yscale,by=yscale/5),
				axisSeq=axisSeq,
				dataStep=as.numeric(step),
				dataSeq=dataSeq,
				yscale=yscale,
				bin_size_in_align=as.numeric(bin_size_align),
				...
		)
		dev.off()
	}
	
	if(pngOutput){
		fileNamePng=paste(prefix,".",fileName,".png",sep="")
		png(fileNamePng)
		heatmap.3(valueTmpGch,
				Rowv=FALSE,
				Colv=FALSE,
				dendrogram= c("none"),
				distfun = dist,
				hclustfun = hclust,
				na.rm=TRUE,
				labRow = "",
				labCol = "",
				#labCol = seq(0-scale,scale,by=step),
				cexCol = 0.05 + 1/log10(dim(valueTmpGch)[2]),
				key=TRUE,
				keysize=1,
				trace="none",
				density.info=c("none"),
				margins=c(10, 8),
				RowSideColors=colBar,
				col=heatMapCols,
				#na.color=par("bg"),
				#col=blue2yellow(75),
				...
		)
		dev.off()
	}
	
	if(dendgramAlone){
		fileNameDend=paste("Dendgram.",prefix,".",fileName,".pdf",sep="")
		pdf(fileNameDend)
		plot(x)
		dev.off()
	}
	returnResult<-list(rowNames=rownames(valueTmpGch), subClusterOrder=subClusterOrder)
	return(returnResult)
	
	
}


generateHeatmapForNOMeSeqInFirstMatrixClusterOrder<-function(fileNames,prefix="heatmap",outputAllInOne=FALSE, sampleLen=101, scale=1000, step=20,rowOrder=NULL, pngOutput=FALSE,pdfOutput=TRUE, ...){
	i=0
	if(outputAllInOne){
		
	}
	else{
		for(i in c(1:length(fileNames))){
			if(i==1){
				rowOrder<-generateHeatmapForNOMeSeq(fileNames[i],rowOrder=NULL, ...)
			}
			else{
				c<-generateHeatmapForNOMeSeq(fileNames[i],rowOrder=rowOrder, ...)
			}
		}
	}
	
	return(rowOrder)
	
}


clara.variance<-function (x, k, metric = "euclidean", stand = FALSE, samples = 5, 
		sampsize = min(n, 40 + 2 * k), trace = 0, medoids.x = TRUE, 
		keep.data = medoids.x, rngR = FALSE, pamLike = FALSE) 
{
	if (inherits(x, "dist")) 
		stop("'x' is a \"dist\" object, but should be a data matrix or frame")
	x <- data.matrix(x)
	if (!is.numeric(x)) 
		stop("x is not a numeric dataframe or matrix.")
	n <- nrow(x)
	if ((k <- as.integer(k)) < 1 || k > n - 1) 
		stop("The number of cluster should be at least 1 and at most n-1.")
	if ((sampsize <- as.integer(sampsize)) < max(2, k + 1)) 
		stop(gettextf("'sampsize' should be at least %d = max(2, 1+ number of clusters)", 
						max(2, k + 1)), domain = NA)
	if (n < sampsize) 
		stop(gettextf("'sampsize' = %d should not be larger than the number of objects, %d", 
						sampsize, n), domain = NA)
	if ((samples <- as.integer(samples)) < 1) 
		stop("'samples' should be at least 1")
	jp <- ncol(x)
	namx <- dimnames(x)[[1]]
	if (medoids.x) 
		ox <- x
	else if (keep.data) 
		stop("when 'medoids.x' is FALSE, 'keep.data' must be too")
	if (stand) 
		x <- scale(x, scale = apply(x, 2, meanabsdev))
	if (keep.data) 
		data <- x
	if ((mdata <- any(inax <- is.na(x)))) {
		jtmd <- as.integer(ifelse(apply(inax, 2, any), -1, 1))
		valmisdat <- 1.1 * max(abs(range(x, na.rm = TRUE)))
		x[inax] <- valmisdat
	}
	else rm(inax)
	doDUP <- nzchar(dup <- Sys.getenv("R_cluster_clara_dup")) && 
			isTRUE(as.logical(dup))
	if ((trace <- as.integer(trace))) 
		cat(sprintf("calling .C(cl_clara, ..., DUP = %s):\n", 
						doDUP))
	res <- .C(cl_clara, n, jp, k, clu = as.double(x), nran = samples, 
			nsam = sampsize, dis = double(1 + (sampsize * (sampsize - 
									1))/2), mdata = as.integer(mdata), valmd = if (mdata) rep(valmisdat, 
								jp) else -1, jtmd = if (mdata) jtmd else integer(1), 
			ndyst = as.integer(if (metric == "manhattan") 2 else 1), 
			as.logical(rngR[1]), as.logical(pamLike[1]), integer(sampsize), 
			integer(sampsize), sample = integer(sampsize), integer(k), 
			imed = integer(k), double(k), double(k), double(k), avdis = double(k), 
			maxdis = double(k), ratdis = double(k), size = integer(k), 
			obj = double(1), avsil = double(k), ttsil = double(1), 
			silinf = matrix(0, sampsize, 4), jstop = integer(1), 
			trace = trace, tmp = double(3 * sampsize), itmp = integer(6 * 
							sampsize), DUP = doDUP)
	if (res$jstop) {
		if (mdata && any(aNA <- apply(inax, 1, all))) {
			i <- which(aNA)
			nNA <- length(i)
			pasteC <- function(...) paste(..., collapse = ",")
			stop(ngettext(nNA, sprintf("Observation %d has *only* NAs --> omit it for clustering", 
									i[1]), paste(if (nNA < 13) 
												sprintf("Observations %s", pasteC(i))
											else sprintf("%d observations (%s ...)", nNA, pasteC(i[1:12])), 
									"\thave *only* NAs --> na.omit() them for clustering!", 
									sep = "\n")), domain = NA)
		}
		if (res$jstop == 1) 
			stop("Each of the random samples contains objects between which\n", 
					" no distance can be computed.")
		if (res$jstop == 2) 
			stop("For each of the ", samples, " samples, at least one object was found which\n could not", 
					" be assigned to a cluster (because of missing values).")
		stop("invalid 'jstop' from .C(cl_clara,.): ", res$jstop)
	}
	res$clu <- as.integer(res$clu[1:n])
	sildim <- res$silinf[, 4]
	disv <- res$dis[-1]
	disv[disv == -1] <- NA
	disv <- disv[upper.to.lower.tri.inds(sampsize)]
	class(disv) <- dissiCl
	attr(disv, "Size") <- sampsize
	attr(disv, "Metric") <- metric
	attr(disv, "Labels") <- namx[res$sample]
	res$med <- if (medoids.x) 
		ox[res$imed, , drop = FALSE]
	if (!is.null(namx)) {
		sildim <- namx[sildim]
		res$sample <- namx[res$sample]
		names(res$clu) <- namx
	}
	r <- list(sample = res$sample, medoids = res$med, i.med = res$imed, 
			clustering = res$clu, objective = res$obj, clusinfo = cbind(size = res$size, 
					max_diss = res$maxdis, av_diss = res$avdis, isolation = res$ratdis), 
			diss = disv, call = match.call())
	if (k > 1) {
		dimnames(res$silinf) <- list(sildim, c("cluster", "neighbor", 
						"sil_width", ""))
		r$silinfo <- list(widths = res$silinf[, -4], clus.avg.widths = res$avsil, 
				avg.width = res$ttsil)
	}
	if (keep.data) 
		r$data <- data
	class(r) <- c("clara", "partition")
	r
}


